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ABSTRACT 

We develop a simple model for estimating the mass growth histories of dark matter 
halos. The model is based on a fit to the formation time distribution, where formation is 
defined as the earliest time that the main branch of the merger tree contains a fraction 
/ of the final mass M. Our analysis exploits the fact that the median formation time as 
a function of / is the same as the median of the main progenitor mass distribution as 
a function of time. When coupled with previous work showing that the concentration 
c of the final halo is related to the formation time tf associated with / ~ 0.04, 
our approach provides a simple algorithm for estimating how the distribution of halo 
concentrations may be expected to depend on mass, redshift and the expansion history 
of the background cosmology. We also show that one can predict log 10 c with a precision 
of about 0.13 and 0.09 dex if only its mass, or both mass and tt are known. And, 
conversely, one can predict log 10 tf from mass or c with a precision of 0.12 and 0.09 dex, 
approximately independent of /. Adding the mass to the c-based estimate does not 
result in further improvement. These latter results may be useful for studies which 
seek to compare the age of the stars in the central galaxy in a halo with the time the 
core was first assembled. 
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1 INTRODUCTION 

One main focus of modern observational cosmology is the 
study of the uabundance and structural properties of galaxy 
clusters: typically, the abundance is expressed as a func- 
tion of mass M, and the structural property of most inter- 
est is the central concentration c of the mass density pro- 
file. There are various ongoing and planned surveys which 
are capable of finding thousands of clusters using well- 
controlled/calibrated selection criteria, so a number of large 
homogeneous samples of galaxy clusters will be available in 
the near future (Laureijs et al. 2011). Such samples have the 
potential to strongly constrain the fundamental parameters 
of a given cosmological model (e.g. ACDM). This is in large 
part because of numerical simulations (Navarro et al. 1997, 
e.g.) which have shown that, at fixed mass, halos are less 
concentrated if the matter density, Q, m> is lower. In addition, 
for a given cosmology, there is a reasonably tight monotonic 
correlation between M and c, although this M — c relation 
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depends on how c is estimated. Recently Prada et al. (2011) 
report that, if c is determined by fitting to the circular ve- 
locity rather than the density profile, then the c—M relation 
is no longer monotonic. Provided one knows which method 
has been used to estimate c, measuring the M — c relation 
at a range of redshifts provides important constraints on the 
background cosmology (Eke et al. 2001; Dolag et al. 2004; 
Neto et al. 2007; Maccio et al. 2008).) The tightness of these 
constraints rests on accurately estimating the mass and con- 
centration of each cluster. 

In all cases, these relations are thought to be a con- 
sequence of the fact that the structure of a halo, and its 
concentration c in particular, depend upon its assembly his- 
tory (Navarro et al. 1996; Bullock et al. 2001; Wechsler et al. 
2002; Zhao et al. 2003b, 2009; Giocoli et al. 2010a). Roughly 
speaking, c is related to the ratio of the background density 
at the time at which the mass within the core of the halo 
was first assembled to that at which the halo was identi- 
fied (e.g., today). Lower formation redshifts typically imply 
smaller concentrations, but the exact dependence is a func- 
tion of the expansion history of the background cosmology. 
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Recently Zhao et al. (2009) have provided a prescrip- 
tion for relating a halo's concentration to its mass assembly 
history. This prescription requires knowledge of the time at 
which a halo had assembled 4% of its mass: 



= 4 



1 + 



ti 

3.75to.o4 



1/8 



(1) 



where tf denotes the time at which the halo had assembled 
a fraction / of its mass (in what follows we will use indis- 
tinctly both t\ and to to denote the present day time, at 
which the halo had assembled all its mass). If this remark- 
able relation were exact, then it would imply that the M — c 
relation, and its dependence on cosmology, comes entirely 
from the fact that to.oi/tx has a distribution which depends 
on halo mass (recall that massive halos assemble their mass 
later, so they typically have larger values of to.oV^i)- This 
mass dependence itself depends on cosmology, as does the 
conversion between time and the observable, redshift. 

Unfortunately, estimating the mass and cosmology de- 
pendence of tf/ti is not straightforward. Analytic arguments 
for this are only available for / > 1/2 (Lacey & Cole 1993; 
Nusser & Sheth 1999), and these are not particularly accu- 
rate (Giocoli et al. 2007). So the main goal of the present 
work is to provide a simple prescription that is accurate even 
when /< 1. 

Although the close relation between c and assembly his- 
tory is one of the main motivations of the present work, we 
note that there is also interest in using the distribution of 
galaxy velocity dispersions (Sheth et al. 2003) and the weak 
lensing signal to constrain cosmological parameters (Fu et al. 
2008; Schrabback et al. 2010). Both these observables are 
expected to be related to halo mass and concentration, so 
these studies will also benefit from a better understanding 
of the formation history of dark matter haloes. 

In addition, there is considerable interest in quantifying 
the difference between the time when the stars in the central 
galaxy in a halo formed and the time when those stars were 
first assembled into a single object. Whereas crude estimates 
of the stellar age can be made from the observed colors, there 
are currently no estimates of the halo assembly time. A by- 
product of our study of the joint distribution of halo mass, 
concentration and assembly history is an estimate of how 
precisely the mass and concentration of a halo predict its 
assembly time. 

This paper is organized as follows. In Section 2 we de- 
scribe the close connection between the formation time dis- 
tribution and the mass of the main progenitor on which our 
model is built. Section 3 describes our simulation-based es- 
timate of the formation time distribution for arbitrary /. 
It also quantifies correlations between halo mass, concentra- 
tion and tf, providing estimates of how precisely the con- 
centration can be estimated from knowledge of the mass and 
assembly history, as well as how well the mass and concen- 
tration of a halo predict its formation time. In Section 4 we 
describe how to use our fitting formula to generate Monte- 
Carlo merger histories from which to estimate halo concen- 
trations. Section 5 quantifies how halo assembly histories 
depend on cosmological parameters, and a final section sum- 
marizes. 



2 FORMATION AND THE MASS OF THE 
MAIN PROGENITOR 

The mass which makes up a halo of mass Mo at zo will be in 
many smaller pieces at z > zo . We will follow common prac- 
tice in calling these pieces the 'progenitors', and we will use 
N(m, z\Mo, zo) to denote the 'mass function' of progenitors: 
i.e., the mean number of progenitors of Mo that, at z > Zq, 
had mass m. We will use rni(z) to denote the most massive 
of the set of progenitors present at z. This 'most massive pro- 
genitor' may be different from what we will call 'the main 
progenitor', raiipfz), which is the most massive progenitor 
of the most massive progenitor of the most massive progeni- 
tor ... and so on back to some high redshift z > z$. Indeed, it 
is almost certainly true that mi(z) ^ jtimp(z), with equality 
only guaranteed when m\(z) > Mo/2. We quantify the dif- 
ferences between these two quantities in the Appendix. But 
because it is jtimp which is relevant to equation (1) we will 
use it in what follows. 

The formation time of a halo, Pf(> Zf|Mo,Zo), is usu- 
ally defined as the earliest time when iriMpfz) > fMo- (Note 
that this is different from defining formation as the earliest 
time that a single progenitor contains a fraction / of the final 
mass, which would correspond to requiring mi(z) > fMo. 
We discuss this choice in the Appendix.) If pmp (m, z\Mo, zo) 
denotes the probability that raiipfz) = m, then 



p f (> z f \Mo,z ) =pmp(> fM ,Zf\M ,Zo), 



(2) 



where the subscript / on the left hand side is a reminder that 
the shape of the distribution depends on /. This is a simple 
generalization, to arbitrary values of /, of the argument in 
Lacey & Cole (1993) who assumed that / = 1/2, and Nusser 
& Sheth (1999) who considered / > 1/2. 

When / ^ 1/2 then ttimp = mi, since a halo cannot 
have more than one progenitor with mass greater than fMo. 
As a result, pMp(m, z\Mo, zo) is simply related to the pro- 
genitor mass function N(m, z\Mo, zq). In particular, 

Pmp{> fM , zj\Mo, zo) = N(> fM , z f \Mo, z ) for / > 1/2. 

(3) 

This is fortunate, since simple analytic estimates of N(> 
fMo, Zf | Mo, zo) are available. Unfortunately, this simplic- 
ity is lost when / < 1/2, and so, although equa- 
tion (2) still holds, there is no analytic estimate for pmp(> 
f Mo, Zf\Mo, zo) in this case. 

In what follows, we will simply measure Pf(> 
Zf\Mo,zo), with Zf determined by requiring m,Mp(zf) ^ 
fMo in simulations, and use it to determine the MP distri- 
bution at any given z. Namely, whereas the usual procedure 
(when / ^ 1/2) is to estimate the formation time distribu- 
tion by differentiating the right hand side of equation (2) 
with respect to Zf, we will instead estimate the MP distri- 
bution (when / < 1/2) by differentiating the left-hand side 
with respect to / (at fixed z). From this, one can derive an 
estimate of the mean value of m,MP at each z. But note that 
this mean value can be got more directly from equation (2) 
by using the fact that 

M 

dmMPPMp(> m M p, z\M , z ) = (mmp{z)\Mo, z ) . 

(4) 

Since the median ?timp at z is easy to estimate - it is that 
/ at which equation (2) equals 1/2 - comparison of the me- 
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dian and the mean provide a simple estimate of how skewed 
PMp(m M p, 2 1 M , z ) is. 

In what follows, we will be more interested in the me- 
dian relation than in the mean. This is because the mean 
formation redshift as a function of the required mass of the 
main progenitor, (zf\f), traces out a different curve than 
does the mean main progenitor mass as a function of red- 
shift (f\z) (Nusser & Sheth 1999). In contrast, the median 
Zf as a function of / traces out the same curve as does the 
median / as a function of z: this convenient property of the 
median relations, which equation (2) makes obvious, has not 
been highlighted before. 

2.1 Scaled units 

In principle, the distribution p/(> Zf\Mo, Zq) depends on /, 
Mo and zo . However, the analytic solution for / ^ 1/2 is 
actually a function of the scaled variable 

= S c (z f )-S c {zo) 
s/SUMo) - 5 (Mo) 

where 5 c (z) is the initial overdensity required for spheri- 
cal collapse at z, extrapolated using linear theory to the 
present time and S(M) is the variance in the linear fluc- 
tuation field when smoothed with a top-hat filter of scale 
R = (3M/47rp) 1//3 , where p is the comoving density of the 
background. When 2/ is rescaled in this way, then the shape 
of p(> Wf\Mo, zo) is essentially the same for all values of Mo 
and zo, although the dependence on / is still strong. E.g., 
for / 1/2, the expected shape of this distribution is 



Pf(w f ) = 




(6) 

(Lacey & Cole 1993; Nusser & Sheth 1999). In what follows, 
we show that in our simulations this remains true for smaller 
values of / as well. Recent work (Paranjape et al. 2011; 
Musso & Sheth 2012) suggests that Wf may not be the ideal 
choice for the scaling variable - this may account for some 
of the small departures from self-similarity that are evident 
in the following plots. This means that we do not need to 
fit a different functional form for each value of /, Mo and 
zo. Rather we need only fit to p(wf), which is a function of 
/ and Wf only. The expression above can become negative 
at / < 1/2, so it is not a viable model for a probability 
distribution in this regime. Therefore, our main goal is to 
provide a fitting formula which works even at small /. 

Although we can use this fitting formula for the full 
distribution p(wf) to estimate the mean (wf), we will be 
more interested in its median wj. This is because the median 
formation redshift for fixed / satisfies 

8 c (z f ) = 5 c {z ) + w f s/S{fM ) - S(M ) (7) 

whereas the median main progenitor mass at fixed z > zo 
satisfies 

S(fMo) = S(Mo) + tM*)-W] 2 (g) 

Wf 

(Because both S(fMo) and w~f depend on /, this last ex- 
pression must be solved to yield /.) In practice, these both 
trace out the same relation, so for the purposes of making 



plots, it suffices to rewrite the expression for the median 
formation redshift as: 

*<=(*/) I S(Mq) I s(fMp) 

JJzV)^ l + Wt iW^)i^{Mo)~ l - (9) 

To an excellent approximation, the left hand side is just the 
ratio of linear theory growth factors at the two redshifts. 
E.g., for an Einstein-de Sitter model (itself an excellent ap- 
proximation for redshifts greater than about 2), this ratio 
would be (do/a/) = (to/tf) 2 ^ 3 . And, for a power-law spec- 
trum the final term on the right hand side is just a function 
of /. Since Wf is also a function of / only, the only mass 
dependence comes from S(M )/8 2 (z ). If we fix the mass, 
then this term decreases as 2o increases, so the distribution 
of (tf/to) is expected to shift towards unity as the redshift 
at which the halos are identified increases. 



3 A FITTING FORMULA FOR THE 
FORMATION TIME DISTRIBUTION 

In what follows, we will use the GIF2 numerical simulation 
(Gao et al. 2004; Giocoli et al. 2008) to calibrate our fitting 
formula for p/(w/). 

3.1 The simulation 

The GIF2 simulation data we analyze below are publicly 
available at http://www.mpa-garching.mpg.de/Virgo. The 
simulation itself is described in some detail by Gao et al. 
(2004). It represents a flat A cold dark matter (ACDM) uni- 
verse with parameters (fi m , o"8, h, Q,bh 2 ) — (0.3, 0.9, 0.7, 
0.0196). The initial density fluctuation spectrum had an in- 
dex n = 1, with transfer function produced by cmbf ast (Sel- 
jak & Zaldarriaga 1996). The simulation followed the evolu- 
tion of 400 3 particles in a periodic cube 110 Mpc//i on a side. 
The individual particle mass is 1.73 x 10 9 Mq/K, allowing us 
to study the formation histories of a large sample of haloes 
spanning a wide mass range. Particle positions and veloci- 
ties were stored at 53 output times, mainly logarithmically 
spaced between 1 + z = 20 and 1 . Halo and subhalo merger 
trees were constructed from these outputs as described by 
Tormen et al. (2004); Giocoli et al. (2007, 2008). 

At each simulation snapshot haloes were identified using 
a spherical overdensity criterium. At each snapshot the den- 
sity of each particle is assumed to be inversely proportional 
to the cube of the distance to the tenth nearest neighbor. 
We take as centre of the first halo the position of the dens- 
est particle. We then grow a sphere of matter around this 
centre, and stop when the mean density within the sphere 
falls below the virial value appropriate for the cosmologi- 
cal model at that redshift; we used the fitting formulae of 
Eke et al. (1996) for the virial density. All the particles in- 
side this sphere are assigned to the new forming halo and 
are removed from the list. Starting from the second densest 
available particle we build the second halo, and so on until 
all the particles have been scanned. Particles not assigned 
to any collapsed objects are defined as field particles. 

For each halo more massive than 10 11 ' 5 Mq/H (i.e. con- 
taining at least ~ 200 particles) in the five snapshots, cor- 
responding to 20 = 0, 0.52, 1.05, 2 and 4.04 we build up 
the merger history tree as follows. Starting from a halo at 
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Figure 1. Dependence of rescaled formation time distribution p(wj) on the required fraction / of the final halo mass that must have 
been assembled. Different symbols show results for different rcdshifts zo at which the hosts were identified in the GIF2 simulation. The 
dashed curve shows equation (6), modified following Giocoli et al. (2007). The solid curve shows our equation (10). 



z — Zo, we define its progenitors at the previous output, 
z = zo + dzi , as all haloes containing at least one particle 
that at z = zo will belong to that halo. The "main progeni- 
tor" at zo + dzi is defined as the one which provides the most 
mass to the halo at Zq. Then we repeat the same procedure, 
now starting with the main progenitor at z — zo + dzi and 
identify its main progenitor at 2 = zo + dzi + dz2 , and so on 
backward in time. Since we always follow the main progeni- 
tor halo, the resulting merger tree consists of a main trunk, 
which traces the main progenitor back in time, rriMp(z), and 
of "satellites", which are all the progenitors which, at any 
time, merge directly onto the main progenitor. 

In what follows, to guarantee a well-defined sample of 
relaxed haloes, we excluded all halos for which mMp(z) ex- 
ceeded the final host halo mass by more than 10%, as well 
as all halos which were, in fact, unbound. 



3.2 The fitting formula 

Figure 1 shows the cumulative distribution of rescaled for- 
mation redshifts. Each panel shows results for a different 
value of /, as indicated in the axis labels / = 0.9,0.5,0.1 
and 0.01. In each panel, symbols show the measurements 
for haloes more massive than 10 1 h~ Mq identified at five 
different redshifts, approximately zo = 0,1/2,1,2, and 4. 
Clearly, scaling Zf to Wf produces curves which are approx- 
imately independent of redshift. 

The solid line shows the result of fitting 

p O-/) = ^f 7 (10) 

e / ' + OLf — 1 

to the measurements. This returns a value a/ for each /, 
which is shown by the symbols in Figure 2. The a/(/) rela- 
tion is well described by 

a f = 0.815 e- 2f3 /f°- 707 , (11) 
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Figure 2. Dependence of ctf, the free parameter in equation (10), 
on the mass fraction / = m^p/Mo required for formation. 



which is shown as the solid curve in each panel. The shaded 
regions around each curve show the 1 and 2a contours of 

Ax 2 , where x 2 («/) = EJ-R - -P(> ™i,f><Xf)% is com- 
puted for each mass and redshift bin. (Equation 10 con- 
tinues to provide a good fit if we define formation using 
mi rather than myip, with the only difference being that 
ctf ~ 0.867 e -2 ^ //°' 8 . This a/ is larger than that for m.MP 
by a factor of 1.06//° : the median Wf is shifted to higher 
redshifts than when formation is defined using the main pro- 
genitor.) 

Figure 1 shows that equation (10) is reasonably accu- 
rate at least between the first and third quartiles of the 
distributions. The associated differential distribution is 



p(w f ) 



It has median 



9P(> w f ) 
dw f 



a/we 



w 2 /2 



" 2 /2 +Q/ - 1] 



Wf = \/2 In (a/ + 1) . 
and the mean, computed via equation (4), is 



( w i) 



1 



a -<*/)* 

k l/2 



(12) 



(13) 



(14) 



For large values of /, oif <C 1, and (wf) — > \Jctf-K 2 /2 
whereas Wf — > (2/n){wf) . In this limit, the distribution is 
rather skewed. However, it becomes less skewed as / in- 
creases. When otf = 1, then (it)/) = \/^/2 — 1.25 whereas 
Wf = y/2 ln(2) = 1.18. And when / = 0.04 then a f = 8 so 
{Wf) = 2.06 whereas Wf = 2.09. 

Since equation (12) is not very accurate in the tails of 
the distribution, we have also fit p(wf) to 



p(w f ) 



a 0.63/ -2 '' 3 --ytwJ 
AoWf e ' f 



(15) 



(i + omr 2/3 )/Pf 

(Wf) = / dw f p{wf)wf 



For this model the mean is 

r(B + l/M 



if 



(16) 



but one must estimate the median numerically. 

Figure 3 compares the differential formation redshift 
distribution for four values of the assembled fraction /, with 
our two fitting functions. The symbols show the measure- 
ments for haloes more massive than 10 11 ' 5 h~ Mq at five dif- 
ferent zo, as in Figure 1. For equation (15) (dashed curves), 
the best-fit parameters 7/ and /3 satisfy 

7/ = 0.12+— i— + 4.3(/-0.24) 2 , (17) 



and 



290 f 1 



P f =3.05e- a67 / + 



3800 



(18) 



where Ao = /3 7^° /F(Bo) is a normalization factor, with 



Figure 4 shows these scalings. Despite the fact that equa- 
tion (15) provides a more accurate description of the mea- 
surements, we have found equation (12) to be more useful 
because it comes with a simple expression for the median 
relation. 



3.3 Test of scaled units 

Figure 5 shows the correlation between the formation red- 
shift (for / = 0.9, 0.5, 0.1 and 0.01) and the host halo mass 
for the same choices of zq as before, in scaled units. I.e., 
the redshift is expressed as S c (zf) — 5 c (zo) and the mass 
as S(fMo) — S(Ado). The different symbols represent differ- 
ent zo (same as Figure 1); they show the median of the 
correlation, i.e., the median of 5 c (zf) — S c (z ) in bins of 
S(fMo) — S(Mo), while the solid lines of the same color 
show the first and third quartiles. The dotted and solid 
lines in each panel show a least-squares fit to the correla- 
tion, and the prediction of our model (equation 13). The 
dashed line in the two top panels shows the prediction, 
iV(> m, Zf\Mo, zo) — 1/2, associated with equation (6). 
Note that, in the format we have chosen to present our re- 
sults, the prediction boils down to a line of slope 1/2, with 
zero-point equal to log 10 w/. 



3.4 The median mass accretion history 

As described above, we have obtained the merger history 
tree for each halo more massive than 10 11 ' 5 /i _1 Mq at zo, for 
four different choices of zo- Figure 6 shows the mass of the 
main halo progenitor, in units of the mass at 20, as a function 
of z. The different panels show results for parent halos of 
different masses, here expressed as fractions of Mq (where 
S(Mo) = <5 2 (0), so Mo corresponds to 8.9 x lO 12 M /ft). 
From top to bottom we show the halo accretion histories 
starting from different redshifts zo- The open circles and 
dashed curve which enclose them represent the median, the 
first and the third quartile of the tree extracted from the 
simulation. The solid curves enclose 95% of the trees. 

Equation (2) shows that, for these median quantities, it 
should not matter if we plot the median Zf at fixed m/M = 
f, or the median / at fixed z. This is indeed the case, and is 
in contrast to what would have happened if we had chosen to 
show the corresponding mean values: the curve which traces 



6 C. Giocoli et al. 2011 




Figure 3. Differential formation time distribution. Open diamonds show the results from the GIF2 simulation for haloes more massive 
than 10 11 ' 5 /i~'Mq identified at five redshifts zq = 0, 0.52, 1.05, 2 and 4.04. The solid and dashed curves show the result of fitting to 
equations (12) and (15). 
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Figure 4. Dependence of the best-fit free parameters in equation (15) on /; curves show the scalings given by equations (17) and (18). 
The bottom panels show the residuals from these relations. 



as a function of / is slightly different from that which 
traces (f\zf) as a function of z. Our equation (2) shows that 
this is a consequence of the mean not being the same as the 
median: the distributions are skewed. 

The dot-dashed curve in each panel shows the predic- 
tion of Zhao et al. (2009): it describes our measurements 
very well. The dashed curve (for zo = only), shows the 
model of van den Bosch (2002). The solid curve in each panel 



shows our model, evaluated as follows: given Mo at redshift 
zo, we compute S(Mo), 6 c (zo), and a table of S(f Mo) for 
any < / < 1. Inverting equation (13) then gives Zf via 
equation (9). This agrees very well with the numerical sim- 
ulation results and therefore also with the model of Zhao 
et al. (2009). But note that our prescription is much easier 
to implement. 
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S[0.1 M vir (z )]-S[M vir (z )] S[0.01 M vir (z )]-S[M vir (zJ] 

Figure 5. Correlation between formation redshift - defined as the first time that the main progenitor contains a fraction / of the final 
mass Mo — and host halo mass, expressed in scaled units. The different panels show results for different choices of / (same as previous 
figure). The symbols, also same as previous figure, show results for different zq. The dotted, dashed, and solid lines show a least-squares to 
the data at all redshifts (open diamonds), the relation associated with equation (6), and the prediction of our simple model (equation 13). 
In this format, both models predict a line of slope 1/2, so they only differ in the value of the zero-point. 



3.5 Correlation between formation times 



So far, we have mainly tested if w/ is indeed a good scaling 
variable for the entire range of /. In what follows, it is in- 
teresting to know if the different formation time definitions 
are correlated. E.g., if a halo is above the median formation 
redshift for one value of /, is it also likely to have formed at 
above average wj for another /? Figure 7 shows that this 
correlation is weak, at least for / = 1/2 and 0.04: When 
wo. 5 — (wo. 5} and u>o.04 — (wo. 04} are normalized by their 
rms values, then the correlation coefficient is ~ 0.25. 



3.6 Halo concentrations from masses and 
formation times 

Equation (1) states that the formation time to. 04 of a halo 
can be converted into an estimate of its concentration. To 
check this, we estimated cw for each halo in GIF2 by fitting 
the matter density profile to the NFW functional form, 

Ps 



p(r\M vir ) = 



r/r s (l + r/r 3 ) 2 



where 



ps = 



M m 



V'lT 



ln(l + c vi 



1 ~\~ Cvir 



(19) 



(20) 



and Cvir = r 3 /R V i r . Figure 8 shows the correlation between 
c V i r and to/to. 04, for various M„; r and zq. The dot-dashed 
curve (same in each panel) shows equation( 1): if cw were 
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Figure 6. Mass growth history of dark matter haloes in the ACDM GIF2 cosmology. The different panels refer to various redshifts and 
host halo masses, as labeled. In each panel, open circles show the median z at fixed m/M, and the curves that enclose them show the 
first and the third quartiles of the distribution defined by the GIF2 sample. These curves are the same as those for the median m/M at 
fixed z: as expected from equation (2). The dot-dashed, dashed and solid curves show the model of Zhao et al. (2009), van den Bosch 
(2002) (valid only for haloes starting from zq = 0) and our equation (9). 
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Table 1. Precision of concentration estimates for halos with 
M vir > 10 13 h- 1 M e . 



-2 



Figure 7. Correlation between formation times, when formation 
is defined as the first time that the main progenitor m^p contains 
0.5 and 0.04 of the total mass. We have scaled the variables by 
their rms values, so the fact that the median ioo.5/°"0.5 a t fixed 
u>o.04/<to.04 is rather different from the median uio.o4/co.04 at 
fixed ioo.5/°"0.5 indicates that the correlation between these two 
is weak. 



completely predicted by to/to. 04 then there should be no 
scatter around this relation. 

In each panel, dots show the halos, and triangles with 
error bars show the median and first and third quartiles in 
log 10 (cw) for narrow bins in to/to. 04- The actual distribu- 
tions of log 10 (c„i r ) and log 10 (£o/fo.04) are shown along the 
left axis and across the bottom. (Note that, for each Mo bin, 
the peak of the distribution of to/to.04 shifts towards unity 
as zq increases, for the reason given following equation 9.) 
Clearly, although equation (1) provides a reasonable descrip- 
tion of the triangles, the scatter around this relation is con- 
siderable. Quantitatively, the scatter in log 10 c V i r is 0.13 dex 
for z = halos with M vir > lO 13 ft~ 1 M . 

The previous subsection showed that to. 5 and to. 04 of a 
halo are only weakly correlated. So it is reasonable to ask if 
some of the scatter in cw can be removed if one includes 
information about to. 5. E.g., one might reasonably expect 
the most concentrated halos to have small values of both 
to. 5 and to. 04. We have found that the typical concentration 
scales as 



4.23 + 



^0.04 / \I0.5 

(21) 

For halos more massive than 10 13 h 1 Mq, the rms scatter 
around this relation (0.1 dex) is indeed reduced with re- 
spect to the case in which only to. 04 is known (0.12 dex). 
Evidently, the remaining scatter is due to other processes. 
Table 1 summarizes how this scatter decreases as more in- 
formation is added. 

One of the main goals of this paper is to provide a simple 
algorithm for generating the M — c relation. Figure 9 shows 
this relation in GIF2. The open circles show the median of 
the distribution at fixed mass, and the error bar encloses 



rms(log 10 


C-vir | ^^vir) 
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rms(log 10 


Cinr|T0.04, Mvir) 


0.120 


rms(log 10 


Cvir\"T0.5, TO. 04, M v i r ) 


0.100 



Table 2. Precision of formation time estimates for halos with 

M vir > 10 13 r 1 ^. 
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rms(-ro.5|<w) 


0.088 
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0.111 
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0.109 
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the first and third quartiles. Different panels show results 
for different redshifts. The dotted line shows a least-squares 
fit to the data, while the dot-dashed and solid curves show 
the relations which come from inserting the mass accretion 
histories of Zhao et al. (2009) and our equation (9) into 
equation (1). Both approaches capture the mass dependence 
of the concentration, and its flattening with redshift, very 
well. But note that our estimate of the value of to.oi/t\ that 
is appropriate for a given mass bin (e.g., our equation 9), is 
much more straightforward to implement. 

The M — c relation for masses of order 1O 6 M0 is needed 
for modeling the milli-lensing effects of galaxy-size halo sub- 
structures (Keeton 2003; Keeton et al. 2003; Sluse et al. 
2003; Metcalf & Madau 2001; Metcalf & Amara 2012). The 
small mass regime is also important for estimating the dark 
matter annihilation signal coming from the integrated ef- 
fect of substructures and field haloes at masses smaller than 
10 6 M Q (Pieri et al. 2005; Giocoli et al. 2008, 2009). How- 
ever, this is below the mass resolution of current numerical 
simulations. Therefore, Figure 10 shows the result of using 
our model (equation (9, where to. 5 an d to. 04 are obtained us- 
ing our mass accretion history model) to estimate the M — c 
relation at six different redshifts down to 3 x W 4 ^ 1 M@. 



3.7 Formation times from masses and 
concentrations 

In real datasets, while we may hope to have a reasonable 
estimate of the mass of the host halo, and perhaps of its 
concentration, we are unlikely to know either to. 04 or to. 5. 
So one might reasonably ask how well tf can be predicted 
given knowledge of either M vir , c V i r or both. 

Since equation (10) gives the distribution of Zf given 
M V i r and zo, it can be manipulated to give the distri- 
bution of t//ti at fixed M V i r - In practice, we can get a 
good idea of the rms scatter around the mean by estimat- 
ing how far on either side of the median value one must 
go before 68% of the total probability has been enclosed. 
These correspond to values of wj = \/2 ln(l + 0.188 a/) and 
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Figure 8. Correlation between host halo concentration and the time at which the system assembled 4% of its mass. The different panels 
(left, central and right) show results for different host halo masses Mo, while the different windows show results for different 20 (same 
as Figure 9). The histograms show the distributions of c v i r and toAo.04> an d the symbols with error bars show the correlation between 
them: the median and first and third quartiles of the distribution of c v i r at each £oAo.04- The dot-dashed curve shows equation (1). 




Figure 9. Median mass concentration relation in the GIF2 simu- 
lation. Open circles show haloes in the simulation at five different 
redshifts Zq, as labeled in the panels. Error bars show the first 
and the third quartiles of the distribution for each mass bin. The 
dotted line shows a least-squares fit to the data; dot-dashed and 
the solid curves show the predictions based on the median mass 
growth history of Zhao ct al. (2009) and our equation (9). In each 
panel we show the rms in logi (c v i r ) at fixed mass, for haloes 
more massive than 10 13 M@/h. 



w f = ln(l + 5.303a/). For / = 0.04 we have a f = 8 so 
10/ = 1.35 and 2.75 respectively. For M = 8.9x lO 12 h _1 M , 
s\M )/5l{zo) = 1, so S c (z f )/S c (z ) = 2.81 and 4.68, making 
(tf/to) = 0.18 and 0.08. Thus, the rms value of log(£//t ) 
should be about 0.16 dex. This is indeed close to the value 
0.17 dex measured directly in the simulations. A similar ana- 
lyis for / = 0.5 returns similarly good agreement. 

The tf — c V i r relation is more complicated. If we 



Figure 10. Mass concentration relation at six different redshifts 
predicted using equation (9). The concentration prediction has 
been extended down to small masses: 3 X 10 4 Mq/H. 



define 77 = log 10 (t//to), then at any Zq, (ro. 5 |c„ ir ) » 
0.175 — 0.503 log 10 Cvir whereas (To.o-ijcw) ~ —0.049 — 
0.825 log 10 c V ir. In both cases, the rms scatter around the 
fit is about 0.1 dex. Adding the mass does not reduce the 
scatter significantly, so we do not report fits to this case. Ta- 
ble 2 summarizes the precision with which a halo's formation 
time can be predicted from its mass and concentration. In 
Figure 11 we show the correlation between the halo con- 
centration and Tf, where / = 0.5 on the left and / = 0.04 
on the right panel, for haloes with M v i r ^ 10 13 Mq//i and 
at any considered zo . The data points with error bars show 
the median, with the quartiles, of both the concentration at 
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Figure 11. Correlation between concentration and Tf (/ = 0.5 on the left and / = 0.04 on the right panel). The data show simulation 
measurements at any zq for all haloes more massive than 10 13 M^/h. Open triangles and circles show the median with the quartiles of 
respectively c v i r at fixed Tf and vice versa. The solid lines represent the least-squares fit to the simulation data points in both axsis. 

fixed tf and vice versa. The solid lines show the leas-squares enclose 95% of the trees. The solid line shows our analytical 
fit to the simulation data in both axsis. prediction (equation 9). 



4 MONTE-CARLO REALIZATIONS OF MASS 
ACCRETION HISTORIES 

We have used the distribution of formation times discussed 
above to generate Monte Carlo realizations of mass accretion 
histories as follows. In reality, wimp changes continuously; 
the amount by which it changes at any given time is a ran- 
dom variable that is different for each halo. Most merger 
tree algorithms approximate this by taking (typically small) 
discrete steps in redshift, and then allow the mass to vary 
stochastically (Nusser & Sheth 1999, e.g.); others allow large 
steps in time (Moreno et al. 2009). Our approach is rather 
different: we predetermine discrete steps in mass, and at 
each step, we draw a redshift from pmp for that mass. 

The choice of mass-step is critical. Two steps in mass 
/iM(i and /2M0 (with fx > f%) would be strongly corre- 
lated if /1 ~ /2, but mass scales separated by A/// ~ 1 
are approximately uncorrelated (Peacock & Heavens 1990; 
Paranjape et al. 2011). Figure 7 shows that, at least for large 
A///, the distributions are indeed reasonably independent. 

We have found that if we set A log 10 M « —0.3 then 
the distributions p(wf 1 ) and p(wf 2 ), when expressed as func- 
tions of redshift, do not overlap much. Since this mass differ- 
ence corresponds approximately to A/// ~ 0.5, we assume 
that we can take independent picks from these redshift dis- 
tributions, with the condition that for the small fraction 
of picks where 2:2 < Z\ occurs we should set 22 = z%. If 
there were no overlap at all, our approach would correspond 
to having independent jumps in redshift, always having the 
correct distribution of redshifts at the predetermined dis- 
crete mass scales. Figure 12 shows the results. The filled 
triangles and the dashed curve that encloses them indicate 
the median and the first and third quartiles, solid curves 



5 APPLICATIONS 

Extended-Press & Schechter (1974) approach gives the pos- 
sibility of estimating the main halo growth as function of the 
growth factor, the mass variance and the overdensity thresh- 
old for collapse. Once fixed the cosmological model these 
quantities can be translated in mass, redshift and time. For 
this reason the halo formation, growth and concentration de- 
pend on the cosmological parameters adopted (Harker et al. 
2006; Power et al. 2011). In the left panel of the Figure 13 we 
show how much the median redshift of haloes at the present 
time, with Zq=Q, differs with respect to a reference model in 
assembling 90%, 50%, 10% and 1% of their mass. We have 
compared as = 0.8 (bottom) and as = 1 (top) with respect 
to <jg = 0.9. as represents the amplitude of the linear power 
spectrum on a scale of 8 Mcp//i, from the figure we notice 
that universes with a lower value of as tend to form haloes 
more recently. The same is valid if we consider the dark 
energy equation of state parameter w. In the right panel of 
Figure 13 we show the relative change in formation redshifts 
between w — —1.1 and w — —0.9 with respect to w = —1. 

Similarly to the previous figure, in Figure 14 we show 
the relative difference in formation redshifts in universes 
with lower and higher matter density parameter with respect 
to the model of the GIF2 simulation (we have assumed that 
Q m + Qa = 1). In the same figure, the right panel quan- 
tify the change in formation redshift due to the change in 
Hubble constant. 

From the figures above we notice that a bigger difference 
in formation redshifts can be mainly seen, for masses of the 
order of galaxy clusters, in universes with a different power 
spectrum normalization as and matter density parameter 
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Figure 12. Mass growth history of a sample of haloes starting from different redshifts zq. The sample as been built sampling a 
discrete number of mass steps with Alog 10 M = —0.3 for the rescaled formation redshit distributions between zq and z > zq. The filled 
triangles and the curves that enclose them indicate the median and the first and the third quartiles. The solid line is our mass growth 
history model (equation 9). 
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Figure 13. Correlation between the formation redshift and the virial host halo mass at zq = 0, for different choices of the mass fraction 
required for formation (/ = 0.9, 0.5 0.1 and 0.01). Left panel shows the relative change in formation redshift (and so the main halo 
growth and the halo abundance) for <rg = 0.8 and erg = 1 relative to erg = 0.9. Right panel shows the dependence on the equation of 
state parameter w. 
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Figure 14. The same as Figure 13, in this case on the left we show the relative difference in formation redshifts for universes with lower 
and higher density parameter with respect to the GIF2 cosmology, the right with different Hubble constant. 



Q m , a change of 0.1 in w and h do not influence too much 
the formation time definitions. 



6 DISCUSSION AND CONCLUSIONS 

If the formation time of a halo is defined as the first time 
that its main progenitor exceeds a fraction / of the final mass 
M, then, for halos of a given M, there will be distribution 
of formation times which depends on /. We argued that the 
median formation time as a function of / should be the same 
as the median of / as a function of time (equation 2), and 
showed that this was indeed the case (Figure 6). 

We also argued that, to a good approximation, the dis- 
tribution of formation times should be a function of the scal- 
ing variable wf defined in equation (5), even when / < 1/2. 
The distribution of scaled formation times p(wf) is indeed 
approximately independent of M and zo (Figure 1). This dis- 
tribution is well-described by equation (10). Its median has 
a simple form. When combined with a relation which relates 
formation time to halo concentration (equation 1), this pro- 
vides a simple way to estimate the concentration of a halo. 
Including more information about the mass accretion his- 
tory (e.g. equation 21) helps predict the concentration more 
precisely (Table 1). 

We also studied how well a halo's formation time could 
be predicted from knowledge of its mass and concentration. 
If the concentration is known, then the mass adds little new 
information about the formation time (Table 2). When the 
formation time is expressed in units of the time ti at which 
the halo was identified, then the concentration estimates 
log 10 (to.s/ii) with rms 0.09 dex, and log 10 (to.o4/ti) with 
rms 0.11 dex. It may be that these scalings will find use in 
studies which compare the mean age of the stars in the cen- 
tral galaxy of a cluster with the time when the mass in the 
central core was first assembled. 

The formation times for sufficiently different / are only 
weakly correlated (Figure 7). We argued why this is not 
unexpected (Section 4) and then used this fact to devise 



a simple Monte-Carlo algorithm for generating mass accre- 
tion histories. The algorithm is rather accurate (Figure 12), 
and allows one to generate realistic mass accretion histo- 
ries over a large range of masses and redshifts. We ex- 
pect it to be useful for studies of how the M — c rela- 
tion depends on the background cosmological model. Al- 
gorithms for estimating the median main halo progenitor 
mass growth history and the concentration-mass relation 
are available here: http://cgiocoli.wordpress.com/research- 
interests / concentrationmah. 
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APPENDIX A: APPENDIX 

The main text defined halo formation using the main pro- 
genitor mMP, rather than the most massive progenitor mi. 
In principle, predicting the distribution of mi(z), the most 
massive of the set of progenitors at z, is just an extreme 
value statistics problem. In practice, calculating this distri- 
bution is complicated by the fact that the partition of Mo 
into mi , • • • , rn n does not involve n independent picks from 
the same distribution. (E.g., for white-noise initial condi- 
tions, the joint distribution of the n picks is known in closed 
form - see Sheth & Lemson (1999) - and it is not simply 
the product of n identical distributions.) Similarly, there is 
no simple expression for the distribution of mMp(z). 
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Figure Al. Distribution of scaled formation times for dark matter haloes identified at zq = 0, when formation is defined as the earliest 
time that jtimp ^ fMg (circles) or as mi ^ fMo (crosses) for different values of /. The solid curve shows our simple fitting function 
(equation 10) with af given by equation (Al) rather than (11). Note that the two definitions coincide for / ^ 1/2, but that otherwise 
mi yields higher formation redshifts. 



Figure Al shows the cumulative formation time distri- 
bution of the sample of haloes identified at zq — 0. The 
open circles (same as Figure 1) show the measurements on 
the merger-tree where the main halo progenitor has been 
defined among the most contributing between two consec- 
utive time steps, while for the crossed data we defined the 
main halo as the most contributing progenitor at zo = 0. 
The dashed and the solid curve are the same as in Figure 1 
and represent our fit and the modified model by Nusser & 
Sheth (1999). As expected, the distributions are essentially 
identical when / = 0.9 and 0.5. However, at smaller values 
of /, the difference between the two definitions increases. 

We can quantify this difference by noting that, if we fit 
equation (10) to this distribution, then the required scaling 
of af with / becomes 



atf = 0.867 e" 2/ ' // 



(Al) 



Since this aj is larger than that in the main text by a fac- 
tor of 1.06// 01 , equation (13) shows that the median Wf is 
shifted to higher redshifts than when formation is defined 
using the main progenitor. 
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